########################
# 14_brazil_pressure.R #
########################

#####################################
## Pressure model, vector notation ##
#####################################


rm(list=ls())

library (runjags)
library (coda)
library (reshape2)
library (random)
library (scales)

print(version)

# platform       x86_64-apple-darwin17.0     
# arch           x86_64                      
# os             darwin17.0                  
# system         x86_64, darwin17.0          
# status                                     
# major          4                           
# minor          2.3                         
# year           2023                        
# month          03                          
# day            15                          
# svn rev        83980                       
# language       R                           
# version.string R version 4.2.3 (2023-03-15)
# nickname       Shortstop Beagle            



# Set working directory
savePathRep <- "/all_data/"
setwd (savePathRep)


RandomSeed<- read.csv("brazil_seeds_1MPE.csv")

##Load Initial values for the mean of ideal points and item parameters
load("brazil_initial_values.Rda")


for (i in 1:11){
	for (j in 2:3){
if(i==11) sn<- "mean" else sn<- i

# Do not change the following seed
randomseed<- RandomSeed[i,j]
set.seed(randomseed)

load (paste0("brazil_rollcall_long_format_", sn, ".Rdata"))

##Construct additional variables
Procedure.PT<- Vote.Features$PT.leader.request.procedure.x * PT.dummy*Pressure$pressure
Procedure.PSDB<- Vote.Features$PSDB.leader.request.procedure.x * PSDB.dummy*Pressure$pressure
Procedure.DEM<- Vote.Features$DEM.leader.request.procedure.x * DEM.dummy*Pressure$pressure
Procedure.PMDB<- Vote.Features$PMDB.leader.request.procedure.x * PMDB.dummy*Pressure$pressure
Procedure.PP<- Vote.Features$PP.leader.request.procedure.x * PP.dummy*Pressure$pressure
Procedure.PR<- Vote.Features$PR.leader.request.procedure.x * PR.dummy*Pressure$pressure
Procedure.leader.pressure<- Procedure.PT+Procedure.PSDB+Procedure.DEM+Procedure.PMDB+Procedure.PP+Procedure.PR
Procedure.leader<- pmax(Vote.Features$PT.leader.request.procedure.x, Vote.Features$PSDB.leader.request.procedure.x, Vote.Features$DEM.leader.request.procedure.x, Vote.Features$PMDB.leader.request.procedure.x, Vote.Features$PP.leader.request.procedure.x, Vote.Features$PR.leader.request.procedure.x)

pcb<- ifelse(Vote.Features$pcb.request.x==1, 1, 0) 
R28 <- PT.dummy*pcb*Pressure$pressure
R29 <- PSDB.dummy*pcb*Pressure$pressure
R30 <- DEM.dummy*pcb*Pressure$pressure
R31 <- PMDB.dummy*pcb*Pressure$pressure
R32 <- PR.dummy*pcb*Pressure$pressure
R33 <- PP.dummy*pcb*Pressure$pressure

 potential.fix.bill<- c(25, 41, 42,71, 240, 241)


###
 
pressure.only.v="model { 
   for (i in 1:n.obs) {
       y[i] ~ dbern(pi[i])
       probit(pi[i]) <- beta[vote[i]]*theta[dep[i]] + eta[vote[i]]*delta[dep[i]]  - alpha[vote[i]] + lambda[1]*P1[obs[i]] +
                        lambda[2]*P2[obs[i]] + lambda[3]*P3[obs[i]] + lambda[4]*P4[obs[i]] +
                        lambda[5]*P5[obs[i]] + lambda[6]*P6[obs[i]] 
                        + lambda[7]*P7[obs[i]]
                        +lambda[8]*P8[obs[i]] + lambda[9]*P9[obs[i]] + lambda[10]*P10[obs[i]] +
                        lambda[11]*P11[obs[i]]+lambda[12]*P12[obs[i]]
                        + lambda[13]*P13[obs[i]]
                        +lambda[14]*P14[obs[i]] + lambda[15]*P15[obs[i]] + lambda[16]*P16[obs[i]] +
                        lambda[17]*P17[obs[i]] +lambda[18]*P18[obs[i]]
                      +lambda[19]*P19[obs[i]]
                      +lambda[20]*P20[obs[i]]
                      +lambda[21]*P21[obs[i]]
                       +lambda[22]*P22[obs[i]]
                       +lambda[23]*P23[obs[i]]
                      +lambda[24]*P24[obs[i]]
                      +lambda[25]*P25[obs[i]]
                      +lambda[26]*P26[obs[i]]

   } 
   # PRIORS
   for(j in 1:26) { lambda[j] ~ dnorm(0,1) } 
# PRIORS
for (j in 1:n.item){ alpha[j] ~ dnorm(0,1) }   
# Beta (discrimination, dimension 1)
for(j in 1:(fix.bill[1]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[1]] ~ dnorm(0,1) 
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[2]]  ~ dnorm(0,1) #
for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[3]] ~ dnorm(0,1)
for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[4]] <- -1.3
for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[5]] <- 1.6
for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[6]] <- -3
for(j in (fix.bill[6]+1):n.item) { beta[j]  ~ dnorm(0,1) }
# Eta (discrimination, dimension 2)
for(j in 1:(fix.bill[1]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[1]] <- -2
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[2]] ~ dnorm(0,1)
for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[3]] ~ dnorm(0,1)
for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[4]]  <- 2 
for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[5]] <- -1.2#
for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[6]] ~ dnorm(0,1) #
for(j in (fix.bill[6]+1):n.item) { eta[j]  ~ dnorm(0,1) }
# ideal points
for(i in 1:n.legs)  { theta[i] ~ dnorm(0,1) }
for(i in 1:n.legs)  { delta[i] ~ dnorm(0,1) }
## Notes
# In this model we only include pressure, but make no attempt to model abstentions
}"
 
 
y1<- molten.roll.call$rc
y2<- y1
for (i in 1:length(y1)){
	if (is.na(y1[i])==FALSE)
		{if(y1[i]==-5) y1[i]<- NA}
	if (is.na(y1[i])==FALSE)
	{if(y1[i]==-7|y1[i]== -9)
	y1[i]<- NA	}
	
}



#molten.roll.call$rc
pressure.data.vector <- dump.format(list(y=y1
                                         , fix.bill=potential.fix.bill
                                         , n.legs=max(molten.party.member$final.leg.no)
                                         , n.item=max(as.numeric(molten.roll.call$vote))
                                         , n.obs=nrow(molten.roll.call)
                                         , obs=seq(1:nrow(molten.roll.call))
                                         , vote=as.numeric(molten.roll.call$vote)
                                         , dep=molten.party.member$final.leg.no
                                     
                                          , P1=R40
                                         , P2=R41
                                         , P3=R42
                                         , P4=R43
                                         , P5=R44
                                         , P6=R45
                                          , P7=R2
                                         , P8=seniority
                                         ,P9=R3
                                         , P10=committee.chair
                                         , P11=R8
                                         , P12=district.income
                                         , P13=R138
                                         , P14=percentile.rank
                                         , P15=R124
                                         , P16=log.district.M
                                        ,P17=R150
                                        ,P18=molten.roll.call$year
                                       ,P19=Procedure.leader.pressure
                                       ,P20=Procedure.leader
                                       ,P21=R4
                                       ,P22=minister
                                      ,P23=Vote.Features$president.position.x*Pressure$pressure
                                       ,P24= Vote.Features$president.position.x
                                       ,P25=pcb*Pressure$pressure
                                       ,P26=pcb
                                       
      
))

pressure.parameters = c("theta", "delta", "alpha", "beta", "eta", "lambda", "deviance")


 set.seed (randomseed)
eta1 = initials$eta + runif(length(initials$eta), -0.2, 0.2)
beta1 = initials$beta + runif(length(initials$beta), -0.2, 0.2)
eta1[potential.fix.bill[4:5]] <- c(NA, NA)
eta1[potential.fix.bill[1]] <- NA
beta1[potential.fix.bill[4]] <- NA
beta1[potential.fix.bill[5:6]] <- c(NA, NA)
pressure.inits.1 <- function() {
   dump.format(
      list(
 theta = initials$theta + runif(length(initials$theta), -0.2, 0.2),
      delta = initials$delta + runif(length(initials$delta), -0.2, 0.2),
      alpha = initials$alpha + runif(length(initials$alpha), -0.2, 0.2),
      eta = eta1,
      beta = beta1
      , lambda = rnorm(26,0,1)
      ,'.RNG.name'="base::Wichmann-Hill"
      ,'.RNG.seed'= randomseed
      ))
}
 
 set.seed(randomseed)
pressure.model <- run.jags(
   model=pressure.only.v,  # add .small for small sample
   monitor=pressure.parameters,
    data=pressure.data.vector,
     inits=list (pressure.inits.1()),
    thin=30, burnin=4000, sample=100,
   check.conv=FALSE, plots=FALSE)
date()

chainsPress <- mcmc.list(list (pressure.model$mcmc[[1]]))

   # Uncomment next line to save new version
   #save (chainsPress, file=paste0("CJRPress Brazil 26 variables 20230506 30t 4000b 100s rn", randomseed, " sample ", sn,".Rda"))
}
}

